function y=elasticity_stderr(beta,Scenario)

%This is a function to calculate the 5% and 95% percentiles of elasticity using the estimates

N=size(beta,2);
S=size(beta,3);
L=size(beta,1);
option=L-size(Scenario,2)+1;

price_legal=Scenario(1,1);
price_street=Scenario(2,1);
dltime=Scenario(3,2);

beta_price=reshape(beta(4,:,:),N,S);
beta_dltime=reshape(beta(5,:,:),N,S);

w=[eye(option-1) Scenario];


prob=zeros(N,option-1);

for s=1:S

    elaslegal=zeros(N,3);
    elasstreet=zeros(N,3);
    elasinternet=zeros(N,3);
    
for id=1:N
    
    elasl=zeros(1,3);
    elass=zeros(1,3);
    elasi=zeros(1,3);
    
    numer=exp(w*beta(:,id,s));
    probi=(numer/(sum(numer)+1));

    for k=1:3
        elasl(k)=beta_price(id,s)*price_legal*probi(1)*((k==1)-probi(k));
        elass(k)=beta_price(id,s)*price_street*probi(2)*((k==2)-probi(k));
        elasi(k)=beta_dltime(id,s)*dltime*probi(3)*((k==3)-probi(k));
    end
    prob(id,:)=probi';
    
    elaslegal(id,:)=elasl;
    elasstreet(id,:)=elass;
    elasinternet(id,:)=elasi;
end

z(s,:,1)=mean(elaslegal./(ones(N,1)*mean(prob)));
z(s,:,2)=mean(elasstreet./(ones(N,1)*mean(prob)));
z(s,:,3)=mean(elasinternet./(ones(N,1)*mean(prob)));
end
size(z(:,:,1),1);
%y=[std(z(:,:,1));std(z(:,:,2));std(z(:,:,3))];
y=[prctile(z(:,1,1),5) prctile(z(:,1,1),95) prctile(z(:,2,1),5) prctile(z(:,2,1),95) prctile(z(:,3,1),5) prctile(z(:,3,1),95);...
    prctile(z(:,1,2),5) prctile(z(:,1,2),95) prctile(z(:,2,2),5) prctile(z(:,2,2),95) prctile(z(:,3,2),5) prctile(z(:,3,2),95);...
    prctile(z(:,1,3),5) prctile(z(:,1,3),95) prctile(z(:,2,3),5) prctile(z(:,2,3),95) prctile(z(:,3,3),5) prctile(z(:,3,3),95)];